Exercice 1 : Manipulate sf objects and associated data.frames

  1. Load in R the iris1 map layer “iris_31.rds” of the french department called Haute-Garonne (numbered 31). Would the sf::st_read function also work? Why ?

CLUES Use the readRDS function.

library(sf)
iris_31 <- readRDS("../data/iris_31.rds")
# iris_31 <- st_read("../data/iris_31.rds")
It would not work because iris_31 is not a shapefile but a file already R formatted. Simply load it with the readRDS function.
  1. Display the basemap of department 31 with plot(iris_31). What do you notice ?
plot(iris_31)

# We notice that R performs 3 graphs: one graph per variable in the sf file.
  1. What is the role of the sf::st_geometry function? What solution do you propose then?
# ?sf::st_geometry
# This function makes it possible to isolate the geography information of the sf object so that we put aside other variables (here CODE_IRIS, P14_POP and INSEE_COM).
plot(st_geometry(iris_31))
  1. In which projection is the map layer? Map it with another projection.

Information

Test the Azimuthal Equidistant projection with “crs=”+proj=aeqd +lat_0=90 +lon_0=0" to see a clear difference and create a layer called “iris_31_aeqd”.

CLUES Use the sf::st_crs function to guess the projection and sf::st_transform to change the projection.

#?st_crs
st_crs(iris_31)
par(mar = c(0,0,0,0), mfrow = c(1,2))
plot(st_geometry(iris_31))
iris_31_aeqd <- st_transform(iris_31, crs="+proj=aeqd +lat_0=90 +lon_0=0")
plot(st_geometry(iris_31_aeqd))

  1. Calculate the distance matrix between the 5 first iris of department 31. Do you get the same distance matrix if you are working on a layer projected in another projection?

Information

Use map layers called “iris_31” and “iris_31_aeqd”
mat <- st_distance(x = iris_31[1:5,], y = iris_31[1:5,])
mat
Units: m
         [,1]      [,2]     [,3]      [,4]     [,5]
[1,]     0.00 51153.341 56509.22 53137.481 13756.98
[2,] 51153.34     0.000 20957.21  3086.011 49828.91
[3,] 56509.22 20957.212     0.00 11077.630 61863.18
[4,] 53137.48  3086.011 11077.63     0.000 54345.46
[5,] 13756.98 49828.910 61863.18 54345.458     0.00
mat_aeqd <- st_distance(x = iris_31_aeqd[1:5,], y = iris_31_aeqd[1:5,])
mat_aeqd
Units: m
         [,1]      [,2]     [,3]      [,4]     [,5]
[1,]     0.00 57206.922 62306.32 59385.205 13798.28
[2,] 57206.92     0.000 20957.48  3204.855 55368.23
[3,] 62306.32 20957.481     0.00 11076.220 66744.65
[4,] 59385.20  3204.855 11076.22     0.000 59843.30
[5,] 13798.28 55368.225 66744.65 59843.302     0.00
identical(mat,mat_aeqd)
[1] FALSE
# The two matrices are different
  1. Using the layer called “iris_31”, create a new aggregated map layer called “com_31” which corresponds to the municipalities of department 31. Also keep in this new layer the information on the population in each municipality.

Information

The map layer called “iris_31” contains the 5 digit codes of municipalities in its variable INSEE_COM and the 2014 population in its column P14_POP.

CLUES Use the classic functions of the dplyr package: select, group_by et summarize. These functions also work with sf objects.

library(dplyr)
com_31 <- iris_31 %>%
  select(INSEE_COM,P14_POP) %>%
  group_by(INSEE_COM) %>% 
  summarize(P14_POP= sum(P14_POP)) %>% 
  st_cast("MULTIPOLYGON")

plot(st_geometry(com_31))

  1. Using the data contained in “sir_31”, add to this layer the number of restaurants per municipality.

Information

The code of each municipality isn’t in the “sir_31” database. To create it, you have to create a variable called INSEE_COM (5 digits) which concatenates the DEPET (2 digits) and COMET (3 digits) variables.
sir_31 <- readRDS("../data/sir_31.rds")

com_31 <- left_join(com_31,
                    sir_31 %>%
  mutate(INSEE_COM=paste0(DEPET,COMET)) %>% 
  group_by(INSEE_COM) %>% 
   summarize(nb_of_rest= n()) %>% 
  st_set_geometry(NULL),
                    by=c("INSEE_COM"="INSEE_COM")
) %>% 
  mutate(nb_of_rest=ifelse(is.na(nb_of_rest),0,nb_of_rest))
  1. Aggregate all the information present in “com_31” at the level of french intercommunalites (called EPCI).

Information

You have to use the database “table_MAUP.rds” to have a match between the municipality code (CODGEO) and intercommunality code (EPCI).
table_MAUP <- readRDS("../data/table_MAUP.rds") %>% 
  select(CODGEO,EPCI)

epci_31  <- com_31 %>%
  left_join(table_MAUP,by=c("INSEE_COM"="CODGEO")) %>% 
  group_by(EPCI) %>% 
  summarize(P14_POP=sum(P14_POP),nb_of_rest= sum(nb_of_rest)) %>% 
  st_cast("MULTIPOLYGON")

plot(st_geometry(epci_31))

Information

We can notice that the number of restaurants is very low or even equals to zero in most municipalities. We will therefore want to create a new map layer that partially uses the zoning of the municipalities and partially the EPCI zoning.
  1. We want to map the number of restaurants at the level of municipalities if its EPCI (intercommunality) contains more than 100 restaurants and at the level of intercommunalities if the intercommunality contains less than 100 restaurants. Create a map layer that meets this need.

CLUES Start by creating the layers EPCI_less100 and COM_more100 corresponding respectively to the EPCIs which contain less than 100 restaurants and to the municipalities which belong to an EPCI containing more than 100 restaurants. Then, use the do.call(rbind, list(EPCI_less100,COM_more100)) statement to merge them.

EPCI_less100 <- epci_31 %>% filter(nb_of_rest < 100) %>% 
  setNames(c("territory","P14_POP","nb_of_rest","geometry"))

COM_more100 <- left_join(com_31,table_MAUP,by=c("INSEE_COM"="CODGEO")) %>% 
  filter(!EPCI%in%EPCI_less100$territory) %>% 
  select(-EPCI) %>% 
  setNames(c("territory","P14_POP","nb_of_rest","geometry"))

par(mar = c(0,0,0,0), mfrow = c(1,2))
plot(st_geometry(epci_31))
plot(st_geometry(EPCI_less100), col="lightblue",add=TRUE)
plot(st_geometry(epci_31))
plot(st_geometry(COM_more100), col="lightblue",add=TRUE)

mix_31 <- do.call(rbind, list(EPCI_less100,COM_more100))

plot(st_geometry(mix_31), col="lightblue")

Municipalities

Intercommunalities

Both

  1. Using the {cartography} package, simply add to each territory of this map a proportional circle layer related to the number of restaurants.

CLUES The propSymbolsLayer function allows you to draw proportional circles.

library(cartography)
plot(st_geometry(mix_31), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
propSymbolsLayer(mix_31, var = "nb_of_rest", inches = 0.2)

Exercise 2 : Map with cartography and ggplot2 packages

  1. Data preparation:

Information

For the creation of bks et cols, use the getBreaks et carto.pal functions of the {cartography} package. For the creation of the “typo” variable, you can use the cut function and apply the parameters digit.lab = 2 and include.lowest = TRUE.
library(sf)
library(cartography)

fra <- st_read("../data/fra.shp", quiet = TRUE)

mix_31 <- readRDS("../data/mix_31.rds")
mix_31 <- mix_31 %>% mutate(nb_rest_10000inhab = 10000 * nb_of_rest / P14_POP) 

bks <- getBreaks(v = mix_31$nb_rest_10000inhab, method = "quantile", nclass = 6)[-c(1:2)]

cols <- carto.pal("orange.pal", 4)

mix_31 <- mix_31 %>% mutate(typo = cut(nb_rest_10000inhab,breaks = bks, dig.lab = 2, include.lowest = TRUE))
  1. With the help of {cartography} package, make the following map which contains in a choropleth layer the variable “nb_rest_10000inhab” and in a proportional circle layer the variable “nb_of_rest”. Do the same thing with the ggplot2 package.

cartography

ggplot2

With cartography :

par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
plot(st_geometry(mix_31), border="azure")
plot(st_geometry(fra), col="ivory", border = "ivory3", xlim = bb[c(1, 3)], 
     ylim = bb[c(2, 4)],add=TRUE)

choroLayer(mix_31, var = "nb_rest_10000inhab", breaks = bks, col = cols, border = "grey80", lwd = 0.4, 
           legend.pos = "topleft", legend.title.txt = "Number of restaurants\nfor 10000 inhabitants", 
           add = TRUE)
propSymbolsLayer(mix_31, var="nb_of_rest", col=viridis::viridis(1,alpha=0.8),border=NA, legend.pos="left", legend.title.txt = "Number of restaurants", add = TRUE)
layoutLayer(title = "Restaurants", sources = "Insee, 2018", author = "Kim & Tim, 2018", 
            theme = "green.pal", col = "darkred", coltitle = "white", postitle = "center",
            frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

With ggplot2 :

library(ggplot2)

map_ggplot <- ggplot() +
  geom_sf(data = fra, colour = "ivory3",
          fill = "ivory") +
  geom_sf(data = mix_31, aes(fill = typo), colour = "grey80") +
  scale_fill_manual(name = "Number of restaurants\nfor 10000 inhabitants",
                    values = cols, na.value = "#303030")+
   geom_sf(data = mix_31 %>%  st_centroid(),
           aes(size= nb_of_rest), color = viridis::viridis(1,alpha=0.8), show.legend = 'point')+
  scale_size(name = "Number of restaurants",
             breaks = c(1, 500, 2000),
             range = c(0,18))+
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(mix_31)[c(1,3)],
           ylim = st_bbox(mix_31)[c(2,4)]
  ) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(
      title = "Restaurants",
      caption = "Insee, 2018\nKim & Tim, 2018"
  )

plot(map_ggplot)

Exercise 3 : Smooth Density analysis of restaurants in Haute-Garonne

Remarques pour Timothée et Kim

Pour Timothée, élements des fonctions fufun et fufun2 à modifier :

  • Leurs noms…
  • Remplacer bks - 1 par bks dans fufun2 et enlever rasdens <- rasdens + 1
  • Créer un paramètre pour la légende pour pouvoir adapter à la carte (500m pour le cours sur Paris et 5km pour la Haute-Garonne pour l’exo)
  • Dans la fonction fufun, ne pas faire de st_intersection mais un st_intersect afin de ne pas avoir tout le rectangle autour de la Haute Garonne mais seulement les carrés de la grille de Haute-Garonne.
  • Mettre des cex = 1 dans les deux fonctions (fufun et fufun2) pour avoir les mêmes légendes ? Ou, si ça fait pas beau dans la partie cours, rajouter en paramètre des fonctions fufun et fufun2.
  • Dans fufun2 : changer la méthode pour enlever les 0 du calcul des quantiles.
Puis Kim : adapter cet exo 3 avec ces nouvelles fonctions.
  1. Load the dataset “sir_31” used previously and map the more than 4,000 restaurants of department 31 with the {mapview} package. Try using different parameters to customize your map.

Information

For example, you can use the map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … parameters of the mapview function.
library(mapview)
library(sf)
library(cartography)
sir_31 <- readRDS("../data/sir_31.rds")

mapview(sir_31, map.types = "OpenStreetMap", col.regions = "#940000", label = paste(sir_31$L2_NORMALISEE, sir_31$NOMEN_LONG, sep = " - "), color = "white", legend = TRUE, layer.name = "Restaurants in SIRENE", homebutton = FALSE, lwd = 0.5) 
  1. Load the layer “iris_31” used previously and apply the fufun2 function below on the SIRENE data and on the iris of the department 31. Use a 2km cellsize. What observation can you make concerning the readability of the map?

Information

fufun2 <- function(feat, adm, cs = 500, method="equal",nclass=12, author="") {
  grid <- st_make_grid(x = adm, cellsize = cs, what = "polygons")
  x <- st_intersects(grid, feat) 
  gg <- st_sf(n = sapply(X = x, FUN = length), grid)
  gg <- st_intersection(gg,adm) #NEW
  plot(adm$geometry, col = NA, border = NA, main="", bg = "#FBEDDA")
  bks <- c(0,unique(round(getBreaks(gg$n[gg$n>0], nclass = nclass, method = method),0)))
  cols <- (mapview::mapviewGetOption("raster.palette"))(length(bks) - 1)
  choroLayer(gg, var = "n", border = NA, add = TRUE, method = method, breaks= bks, col= cols, legend.pos="topright", legend.values.cex=0.7, legend.title.cex=0.7, legend.title.txt = "resto per ha", legend.nodata=FALSE,legend.values.rnd=0)
  layoutLayer(title = "Count points in a regular grid", scale = 0.5, tabtitle = TRUE, frame = FALSE, sources = "Insee, SIRENE 2018", author = author)
}
iris_31 <- readRDS("../data/iris_31.rds")
fufun2(sir_31, iris_31, 2000)

# La carte est peu lisible en raison du choix de discrétisation
  1. Try to make a more readable map.

CLUES Test the method = “quantile” parameter and vary the number of classes until you get something nice.

fufun2(sir_31, iris_31, 2000,method="quantile",nclass=20)

  1. In order to make the result even more meaningful, we will now make a smoothed map. You can use the fufun function below. Use a 2km sigma and a 2km resolution and the quantile method (n = 12).

Information

library(spatstat)
library(raster)
library(maptools)
fufun <- function(feat, adm, sigma = 100, res = 50, method="equal",nclass=12, author="") {
  bb <- as(feat, "Spatial")
  bbowin <- as.owin(as(adm, "Spatial"))
  pts <- coordinates(bb)
  p <- ppp(pts[, 1], pts[, 2], window = bbowin)
  ds <- density.ppp(p, sigma = sigma, eps = res)
  rasdens <- raster(ds) * 10000 * 10000
  proj4string(rasdens) <- "+init=epsg:2154"
  bks <- unique(round(getBreaks(values(rasdens), nclass = nclass, method = method),0)) 
  cols <- (mapview::mapviewGetOption("raster.palette"))(length(bks) - 1)
  plot(adm$geometry, col = NA, border = NA, main = "", bg = "#FBEDDA")
  plot(rasdens, breaks = bks, col = cols, add = T, legend = F)
  legendChoro(pos = "topright", cex = 0.7, title.cex = 0.7, title.txt = "resto per ha", 
              breaks = bks, nodata = FALSE, values.rnd = 0, col = cols)
  layoutLayer(title = "Smooth density Analysis", scale = 0.5, tabtitle = TRUE, frame = FALSE, sources = "Insee, SIRENE 2018", author = author)
 }
fufun(sir_31, iris_31,  2000, 2000,method="quantile",nclass=12)
  1. Try to improve the map resolution.
fufun(sir_31, iris_31,  2000, 500,method="quantile",nclass=12)

  1. Would we have obtained the same smoothing for this department if we had kept the restaurants of the bordering departments?
No. The borders would have been different.

Exercise 4 : Cartogram

Make a cartogram of department 31 at the level of intercommunalities in proportion to the number of restaurants (SIRENE data).

Information

Load the layer “epci_31.rds” from previous exercises and the {cartogram} package.
epci_31 <- readRDS("../data/epci_31.rds")

library(cartogram)
library(sf)

cartogramme <- cartogram_cont(epci_31, "nb_of_rest", itermax = 5, maxSizeError = 1)
plot(st_geometry(cartogramme), col="lightblue")

Normal borders

Cartogram

Exercise 5 : Linemap

Make a linemap of department 31 at the level of municipalities in proportion to the number of restaurants (SIRENE data).

Information

Load the map layer “com_31.rds” from previous exercises and the {linemap} package.

CLUES Use the two getgrid andlinemap functions of this package. The following parameters work: cellsize = 1750, k = 400 and threshold = 0.01.

com_31 <- readRDS("../data/com_31.rds")

library(linemap)
grid <- getgrid(x = com_31, cellsize = 1750, var = "nb_of_rest")
plot(sf::st_geometry(com_31), col="ivory1", border = NA)
opar <- par(mar=c(0,0,0,0), bg = "ivory2")
linemap(x = grid, var = "nb_of_rest", k = 400, threshold = 0.01,
        col = "ivory1", border = "ivory4", lwd = 0.6, add = TRUE)



reproducibility

sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
[3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
[5] LC_TIME=French_France.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] linemap_0.1.0       cartogram_0.1.0     maptools_0.9-2      raster_2.6-7       
 [5] sp_1.2-5            spatstat_1.56-1     rpart_4.1-10        nlme_3.1-128       
 [9] spatstat.data_1.3-1 mapview_2.3.0       cartography_2.1.1   bindrcpp_0.2       
[13] sf_0.6-3            knitr_1.20          stringr_1.3.1       dplyr_0.7.4        
[17] purrr_0.2.4         readr_1.1.1         tidyr_0.7.2         tibble_1.4.2       
[21] ggplot2_3.0.0       tidyverse_1.1.1    

loaded via a namespace (and not attached):
 [1] satellite_1.0.1      lubridate_1.7.1      webshot_0.5.0        RColorBrewer_1.1-2  
 [5] httr_1.3.1           rprojroot_1.2        tools_3.3.2          backports_1.1.1     
 [9] rgdal_1.2-11         R6_2.2.2             mgcv_1.8-15          rgeos_0.3-24        
[13] DBI_1.0.0            lazyeval_0.2.1       colorspace_1.3-2     withr_2.1.2         
[17] gridExtra_2.3        mnormt_1.5-5         leaflet_1.1.0        rvest_0.3.2         
[21] xml2_1.1.1           scales_0.5.0.9000    classInt_0.1-24      psych_1.7.5         
[25] goftest_1.1-1        digest_0.6.15        foreign_0.8-67       spatstat.utils_1.9-0
[29] rmarkdown_1.8        R.utils_2.5.0        base64enc_0.1-3      pkgconfig_2.0.1     
[33] htmltools_0.3.6      htmlwidgets_1.2      rlang_0.2.1          readxl_1.1.0.9000   
[37] shiny_1.0.5          bindr_0.1            jsonlite_1.5         crosstalk_1.0.0     
[41] R.oo_1.21.0          magrittr_1.5         Matrix_1.2-7.1       Rcpp_0.12.16        
[45] munsell_0.4.3        abind_1.4-5          viridis_0.4.0        R.methodsS3_1.7.1   
[49] stringi_1.1.7        yaml_2.1.18          plyr_1.8.4           grid_3.3.2          
[53] parallel_3.3.2       forcats_0.2.0        deldir_0.1-14        lattice_0.20-34     
[57] haven_1.1.1          tensor_1.5           hms_0.4.0            pillar_1.1.0        
[61] gdalUtils_2.0.1.7    reshape2_1.4.3       codetools_0.2-15     stats4_3.3.2        
[65] unilur_0.3.0.9000    glue_1.2.0           evaluate_0.10.1      modelr_0.1.1        
[69] png_0.1-7            httpuv_1.3.5         foreach_1.4.3        cellranger_1.1.0    
[73] gtable_0.2.0         polyclip_1.9-1       assertthat_0.2.0     mime_0.5            
[77] xtable_1.8-2         broom_0.4.2          e1071_1.6-8          class_7.3-14        
[81] viridisLite_0.3.0    iterators_1.0.8      units_0.6-0         

  1. In French, IRIS is an acronym of ‘aggregated units for statistical information’. Their target sizes are 2000 residents per basic unit.